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1 Introduction 



The dark matter problem remains one of the most puzzhng questions in cosmology. Cos- 
mological analyses reveal that dark matter may be composed of non-baryonic particles, but 
their nature is still to be discovered. Many new physics models try to provide a natural 
solution to the dark matter problem. Supersymmetry in particular offers a stable particle, 
the lightest super symmetric particle (LSP), if i?-parity is conserved, which could be the 
main component of the dark matter in the Universe. The current density of the LSP can 
be calculated and is referred as relic density. Compared to the latest precise WMAP mea- 
surements of the dark matter density [T], relic density can impose stringent constraints on 
the supersymmetric parameters. 

Superlso Relic is an extension of the Superlso program to the calculation of the relic 
density. The program calculates the relic density as well as the flavor physics observables 
using a SUSY Les Houches Accord file (SLHAl [2] or SLHA2 [3]) as input, either generated 
automatically via a call to SOFTSUSY [Ij or ISAJET [5], or provided by the user. The 
calculation can be performed automatically for different supersymmetry breaking scenarios, 
such as the minimal Supergravity Grand Unification (mSUGRA) also called Constrained 
MSSM (CMSSM), the Non-Universal Higgs Mass model (NUHM), the Anomaly Medi- 
ated Supersymmetry Breaking scenario (AMSB) and the Gauge Mediated Supersymmetry 
Breaking scenario (GMSB). 

One of the most important features of Superlso Relic in comparison to the other public 
relic density calculation codes, DarkSusy |6j, IsaRed [7j and Micromegas [3, is that it pro- 
vides the possibility to alter the underlying cosmological model, by modifying for example 
the radiation equation-of-state, the expansion rate or the thermal properties of the Universe 
in the period before Big-Bang nucleosynthesis (BBN), which is experimentally inaccessible 
and remains theoretically obscure. In [OjlTO], it was shown that a modification of the ex- 
pansion rate or of the entropy content of the Universe before BBN can strongly modify the 
calculated relic density and therefore change the relic density constraints on supersymmetric 
parameter space. Superlso Relic makes it possible to evaluate the uncertainties on the 
relic density due to the cosmological model, and inversely, to make prediction on the early 
Universe properties using the particle physics constraints. 

In the following, first the content of the Superlso Relic package will be presented, as 
well as the list of the main routines used for the relic density calculation. The procedure 
to use Superlso Relic will be then explained, and the inputs and outputs of the program 
will be introduced. Finally, some examples of results obtained with Superlso Relic will 
be given. In the Appendices, a description of the formulas and models used for computing 
the relic density will be detailed. 

2 Content of the Superlso Relic package 

Superlso Relic is a mixed C / Fortran program devoted to the calculation of the relic 
density and able to compute many flavor observables. Seven main programs are provided in 
the package, but the users are also invited to write their own main programs. In particular 
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slha.c can scan files written following the SUSY Les Houches Accord formats, and cal- 
culates the implemented observables. The main programs msugra.c, amsb.c, gmsb.c, and 
nuhm.c have to be linked to the SOFTSUSY [1] and/or the ISASUGRA/ISAJET packages, 
in order to compute supersymmetric mass spectra and couplings within respectively the 
mSUGRA, AMSB, GMSB or NUHM scenarios. 

The main steps to compute the observables in Superlso Relic are the following: 

• Generation of a SLHA file with ISAJET or SOFTSUSY (or supply of a SLHA file by the 
user), 

• Scan of the SLHA file, 

• Calculation of the widths of the Higgs bosons with FeynHiggs, 

• Calculation of the squared amplitudes of the annihilation diagrams involved in the 
relic density calculation, 

• Computation of the thermally averaged total annihilation cross section, 

• Solving of the Boltzmann equations and computation of the relic density, 

• Calculation of the flavor physics observables. 

It should be noted that the relic density calculation is performed even if the LSP is a charged 
particle. A theoretical description of the calculation of the thermally averaged total anni- 
hilation cross section can be found in Appendix [X] and the detail of the calculation of the 
relic density in the cosmological standard model is given in Appendix iBl We refer to [llj 
for a complete description of the calculation of the flavor observables. 

The processes involved in the relic density calculation are all the annihilation and co- 
annihilation processes of the type 

i + j^k + l (1) 

where i,j are supersymmetric particles and k,l are Standard Model particles. The number 
of involved processes is more than 3000, and the number of diagrams is even larger. To 
generate all the squared amplitudes, we have written a Mathematica pj] script which calls 
FeynArts [13] and FormCalc [13] and generates the necessary routines for the numerical 
computation of the amplitudes. These routines are part of Superlso Relic and can be 
found in src/relic. They rely on the LoopTools [14] package to initialize FormCalc inter- 
nal routines, as well as FeynHiggs [15] to calculate the widths of the Higgs bosons which is 
performed at two-loop level. 

LoopTools 2 . 4 and FeynHiggs 2.6.5 are included in the Superlso Relic v2 . 5 package, 
and they can be found in src/ contrib. Therefore the user does not need to download these 
programs separately. 

The compilation process of all the needed routines is very long (a few hours), and their 
calculation can take time. Fortunately, all the squared amplitude routines are not neces- 
sary at the same time, as some processes have only negligible effects. Therefore, all the 
squared amplitudes are not computed for a SUSY parameter space point, and a selection is 
performed to save time, as described in Appendix El 
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2.1 Parameter structures 



The package Superlso Relic relies on the definition of a main structure in src/ include .h, 
which is defined as fohows: 

typedef struct parameters 

/* structure containing all the scanned parameters from the SLHA file */ 
{ 

int model; /* mSUGRA = 1, GMSB = 2, AMSB = 3 */ 
int generator; /* ISA JET = 1, SOFTSUSY = 2 */ 

double Q; /* Qmax ; default = M_EWSB = sqrt (m_stopl*mstop2) */ 

double m0,ml2,tan_beta,sign_mu,A0; /* mSUGRA parameters */ 

double Lambda, Mmess,N5,cgrav,m32; /* AMSB, GMSB parameters */ 

double mass_Z,mass_W,mass_b,mass_top_pole,mass_tau_pole; /* SM parameters */ 

double inv_alpha_em, alphas_MZ , alpha, Gfermi ,GAUGE_Q ; /* SM parameters */ 

double charg_Umix [3] [3] , charg_Vmix [3] [3] ,stop_mix[3] [3] ,sbot_mix[3] [3] , 

stau_mix [3] [3] ,neut_mix [6] [6] ,mass_neut [6] ; /* mass mixing matrices */ 

double Min , Ml_Min , M2_Min , M3_Min , At _Min , Ab_Min , Atau_Min , M2Hl_Min , M2H2_Min , 

mu_Min,M2A_Min,tb_Min,mA_Min; /* optional input parameters at scale Min */ 

double MeL_Min , MmuL_Min , MtauL_Min , MeR_Min , MmuR_Min , 

MtauR_Min; /* optional input parameters at scale Min */ 

double MqLl_Min , MqL2_Min , MqL3_Min , MuR_Min , McR_Min , MtR_Min , 

MdR_Min,MsR_Min,MbR_Min; /* optional input parameters at scale Min */ 

double N51 ,N52,N53,M2H1_Q,M2H2_Q; /* optional input parameters (N51...3: GMSB) */ 

double mass_d,mass_u,mass_s ,mass_c ,mass_t ,mass_e ,mass_nue ,mass_mu, 

mass_num,mass_tau,mass_nut ; /* SM masses */ 

double mass_gluon,mass_photon,mass_ZO ; /* SM masses */ 

double mass_hO ,mass_HO ,mass_AO ,mass_H ,mass_dnl ,mass_upl ,mass_stl ,mass_chl , 
mass_bl ,mass_tl ; /* Higgs & superparticle masses */ 

double mass_el ,mass_nuel ,mass_mul ,mass_numl ,mass_taul ,mass_nutl ,mass_gluino , 
mass_chal ,mass_cha2 ; /* superparticle masses */ 

double mass_dnr ,mass_upr ,mass_str ,mass_chr ,mass_b2,mass_t2,mass_er ,mass_mur , 

mass_tau2; /* superparticle masses */ 

double mass_nuer ,mass_numr ,mass_nutr ,mass_graviton, 

mass_gravitino ; /* superparticle masses */ 

double gp,g2,g3,YU_Q,yut [4] ,YD_Q,yub[4] ,YE_Q,yutau[4] ; /* Yukawa couplings */ 
double HMIX_Q , mu_Q , tanb_GUT , Higgs_VEV , mA2_Q , MSOFT_Q , M1_Q , M2_Q , 
M3_Q; /* parameters at scale Q */ 

double MeL_Q , MmuL_Q , MtauL_q , MeR_Q , MmuR_Q , MtauR_Q , MqLl_Q , MqL2_Q , MqL3_Q , MuR_Q , McR_Q , 
MtR_Q,MdR_Q,MsR_Q,MbR_Q; /* masses at scale Q */ 

double AU_Q,A_u,A_c,A_t,AD_Q,A_d,A_s,A_b,AE_Q,A_e,A_mu,A_tau; /* trilinear couplings */ 
/* SLHA2 */ 

int NMSSM , Rpar ity , CPviolat ion , Flavor ; 

double mass_nutau2 ,mass_e2 ,mass_nue2 ,mass_mu2 ,mass_numu2 ,mass_d2 ,mass_u2 , 
mass_s2 ,mass_c2 ; 
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double CKM_lambda , CKM_A , CKM_rho , CKM_eta ; 

double PMNS_thetal2 , PMNS_theta23 , PMNS_thetal3 , PMNS_deltal3 , PMNS_alphal , PMNS_alpha2 ; 
double lambdaNMSSM_Min , kappaNMSSM_Min , AlambdaNMSSM_Min , AkappaNMSSM_Min , 
lambdaSNMSSM_Min , xiFNMSSM_Min , xiSNMSSM_Min ,mupNMSSM_Min , mSp2NMSSM_Min , mS2NMSSM_Min , 
mass_H03 , mass_A02 , NMSSMRUN_Q , lambdaNMSSM , kappaNMSSM , AlambdaNMSSM , AkappaNMSSM , 
lambdaSNMSSM,xiFNMSSM,xiSNMSSM,niupNMSSM,mSp2NMSSM,mS2NMSSM; /* NMSSM parameters */ 
double PMNSU_Q , CKM_Q , MSE2_Q , MSU2_Q , MSD2_Q , MSL2_Q , MSQ2_Q , TU_Q , TD_Q , TE_Q ; 

double CKM[4] [4] ; /* CKM matrix */ 

double H0_mix[4] [4] ,A0_mix[4] [4] ; /* Higgs mixing matrices */ 

double sU_mix[7] [7] ,sD_mix[7] [7] ,sE_mix[7] [7] , sNU_mix[4] [4] ; /* mixing matrices */ 

double sCKM_msq2[4] [4] , sCKM_msl2 [4] [4] , sCKM_msd2 [4] [4] , sCKM_msu2 [4] [4] , 

sCKM_mse2 [4] [4] ; /* super CKM matrices */ 

double PMNS_U[4] [4] ; /* PMNS mixing matrices */ 

double TU[4] [4] ,TD[4] [4] ,TE[4] [4] ; /* trilinear couplings */ 

/* non-SLHA*/ 

double mass_b_lS ,mass_b_pole ,mtmt ; 
double Lambdas ; /* Lambda QCD */ 

/* Flavor constants */ 

double f _B,f _Bs,f _Ds,m_B,m_Bs,m_Ds,m_K,m_Kstar,m_D,life_B,lif e_Bs,life_Ds; 
/* Decay widths */ 

double width_hO , width_HO , width_AO , width_H ; 
} 

parameters ; 

This structure contains all the important parameters and is called by most of the main 
functions in the program. An additional structure specific to the relic density calculation 
is also defined: 

typedef struct relicparam 

/* structure containing the cosmological model parameters */ 
{ 

double ddO,ndd; 
double sdO,nsd; 
double table_ef f [276] [3] ; 

} 

relicparam; 

This structure is used to define the cosmological model based on which the relic density 
calculation is performed. 

2.2 Main routines 

We now review the main routines of the code needed for the relic density calculation. For 
the main procedures related to the flavor observable calculations we refer the reader to [H] . 
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The most relevant C routines are the fohowing: 

• void Init_parain( struct parameters* param) 

This function initiahzes the param structure, setting ah the parameters to 0, apart 
from the SM masses and the value of the strong coupling constant at the Z-boson 
mass, which receive the values given in the PDG08 [IB] . 

• int Les_Houches_Reader (char name [] , struct parameters* param) 

This routine reads the SLHA file named name, and put all the read parameters in 
the structure param. It should be noted that a negative value for param->model in- 
dicates a problem in reading the SLHA file, or a model not yet included in Superlso 
(such as i2-parity breaking models). In this case, Les_Houches_Reader returns 0, 
otherwise 1. 

• int test_slha(char name[]) 

This routine checks if the SLHA file is valid, and if so returns 1. If not, -1 means 
that in the SLHA generator the computation did not succeed {e.g. because of tachy- 
onic particles), -2 means that the considered model is not currently implemented in 
Superlso, and -3 indicates that the provided file is either not in the SLHA format, 
or some important elements are missing. 

• int sof tsusy_sugra(double mO, double ml2, double tanb, double AO, 
double sgnmu, double mtop, double mbot, double alphas_mz, char name[]) 

• int isajet_sugra(double mO, double ml2, double tanb, double AO, 
double sgnmu, double mtop, char name [] ) 

• int sof tsusy_gmsb (double Lambda, double Mmess, double tanb, int N5, 

double cGrav, double sgnmu, double mtop, double mbot, double alphas_mz, 
char name [] ) 

• int sof tsusy_amsb (double mO, double m32, double tanb, double sgnmu, 
double mtop, double mbot, double alphas_mz, char name [] ) 

• int sof tsusy_nuhm(double mO, double ml2, double tanb, double AO, double mu, 
double mA, double mtop, double mbot, double alphas_mz, char name [] ) 

The above routines call SOFTSUSY or ISAJET to compute the mass spectrum cor- 
responding to the input parameters (more details are given in the next sections), and 
return a SLHA file whose name has to be specified in the string name. 
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• void Modelini (struct parameters* param, double relicmass, double maxenergy) 

This routine is an interface between the C routines and the Fortran routines and 
it defines ah the Fortran variables using the C variables. 

• double findrelicmass (struct parameters* param, int *scalar) 

This function determines the LSP mass, and checks if the LSP is scalar (*scalar=l) 
or fermionic (*scalar=0). 

• int Weff (double* res, double sqrtS, struct parameters* param, double relicmass) 

This function calls the Fortran routines and returns the effective annihilation rate 
PVeflF at a given center of mass energy sqrtS, following the procedure described in 
Appendix [Al 

• int Init_relic (double Wefftab[NMAX] [2] , int *nlines_Wef f , struct parameters* param) 

This routine computes for different values of y/s the effective annihilation rates Wes 
needed for the calculation of {av) using the Weff function, and collects them in table 
Wefftab. 

• double sigmav(double T, double relicmass, double Wefftab [NMAX] [2] , int nlines, 
struct parameters* param) 

This function computes the averaged annihilation cross section (av) using the ef- 
fective annihilation rates Wcs collected in table Wefftab. 

• double heff (double T, struct relicparam* paramrelic) 
double sgstar (double T, struct relicparam* paramrelic) 
double geff (double T, struct relicparam* paramrelic) 

These three functions compute respectively h^s, ^fgl and ^efr at the temperature T. 

• double Yeq(double T, struct parameters* param, struct relicparam* paramrelic) 
double dYeq_dT (double T, struct parameters* param, struct relicparam* paramrelic) 

The first function computes Ygg at a temperature T, and the second one its derivative. 

• double Tfo(double Weff tab [NMAX] [2] , int nlines_Weff, double relicmass, 
struct parameters* param, double d, struct relicparam* paramrelic) 
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This function computes the freeze-out temperature following Eq. (|2ip. using the 
Wef f tab generated previously. 

• double relic_density (double Wefftab[NMAX] [2] , int nlines_Weff, 
struct parameters* param, struct relicparam* paramrelic) 
double relic_calculator (char name[]) 

This main procedure computes the relic density using the Wef f tab generated pre- 
viously. relic_calculator is a container function which scans the SLHA file and 
computes the relic density. 

• void Init_cosnioniodel (struct relicparam* paramrelic) 

void Init_modelef f (int model_eff, struct relicparam* paramrelic) 

void Init_dark_density (double ddO, double ndd, struct relicparam* paramrelic) 

void Init_dark_entropy (double sdO, double nsd, struct relicparam* paramrelic) 

These procedures define the cosmological model based on which the relic density is 
computed. Init_cosmomodel has to be called first to initialize the paramrelic struc- 
ture. To alter the QCD equation-of-state as in Appendix [Cl Init_modelef f must be 
called while specifying the model: model_ef f = 1 • • • 5 corresponds respectively to the 
models A, B, B2, B3 and C developed in [T7], and model_eff= to the old model 
formerly used in Micromegas and DarkSusy. If not specified, the model is set by 
default to B (model_ef f = 2). Init_dark_density adds a dark energy density as in 
Eq. (I3ip . with ddO=Kp and ndd=np, and Init_dark_entropy adds a dark entropy 
density as in Eq. (j32p . with sdO=Ks and nsd=n<(. If these routines are not called, no 
additional density will be added, and the calculation will be performed in the standard 
cosmological model. 

• double dark_density (double T, struct relicparam* paramrelic) 
double dark_entropy (double T, struct relicparam* paramrelic) 

double dark_entropy_derivative (double T, struct relicparam* paramrelic) 
double dark_entropy_Gd (double T, struct relicparam* relicparam) 

These functions compute energy and entropy densities needed for the alternative cos- 
mological models described in Appendix iDl 

• int FeynHiggs (char name [] , struct parameters* param) 

This routine calls FeynHiggs to compute the widths and masses of the Higgs bosons 
corresponding to the SLHA file name at the two-loop level, and puts these variables 
in the param structure. 
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The complete list of C procedures implemented in Superlso Relic is available in src/ include .h. 

The Fortran routines can be found in src/relic. They have been generated automati- 
cally by a Mathematica/FormCalc script and they perform the computation of all squared 
amplitudes. Because of the large number of these routines they will not be described further 
here. For the FormCalc specific routines, we refer the reader to the FormCalc manual |14j . 

3 Compilation and installation instructions 

The Superlso Relic package can be downloaded from: 



http : / / superlso . in2p3 . f r/relic 



and is available in two versions: 

• the shared library based version, Superlso Relic Shared, compiles the squared am- 
plitude procedures on-the-fly, if they are needed. The compilation of this version is 
fast, but its running is slow. It is mostly intended for calculation of test-points. 

• the static library based version, Superlso Relic. Here all the squared amplitude 
procedures need to be compiled before running, and therefore the compilation process 
can take a couple of hours. This version computes much faster than the shared library 
version, and is therefore intended for the calculation of a large number of SUSY points. 

The following main directory is created after unpacking: 
superiso_vX . X 

This directory contains the src/ directory, in which all the source files can be found. The 
main directory contains also a Makefile, a README, seven sample main programs (msugra. c, 
amsb . c, gmsb . c, nuhm. c, slha. c, test_modelef f . c and test_standmod. c) and one exam- 
ple of input file in the SUSY Les Houches Accord format (example . lha). The compilation 
options should be defined in the Makefile, as well as the path to the ISAJET isasugra.x 
and SOFTSUSY softpoint.x executable files, when needed. The important compilation 
options to be set are the C compiler name (by default gcc) and the Fortran compiler name 
(by default gf ortran), and their specific flags if needed. Instructions are given in Makefile 
for the Intel compilers and for the Fortran compiler g77. For the Superlso Relic Shared 
version the LD linker name (by default Id) and its specific flags must also be set. 
Superlso Relic is written for a C compiler respecting the C99 standard and a Fortran 
compiler. In particular, it has been tested successfully with the GNU C and GNU Fortran 
Compilers and the Intel C and Intel Fortran Compilers on Linux and Mac 32-bits or 64-bits 
machines, and with the latest versions of SOFTSUSY and ISAJET. Additional information can 
be found in the README flle. 
To compile the library, type 

make 
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This creates libisospin. a in src/ and librelic . a in src/relic, and compiles FeynHiggs 
and LoopTools. Then, to compile one of the seven programs provided in the main directory, 
type 

make name or make name . c 

where name can be msugra, amsb, gmsb, slha, nuhm, test_modelef f or test_standmod. 
This generates an executable program with the . x extension. Note that slha, test_modelef f 
and test_standmod do not need ISAJET or SOFTSUSY programs, but msugra, amsb, gmsb 
and nuhm need at least one of them. 

slha.x calculates the implemented observables, using the parameters contained in the 
SLHA file whose name has to be passed as input parameter. 

amsb.x, gmsb.x, msugra. x and nuhm.x compute the observables, starting first by calcu- 
lating the mass spectrum and couplings thanks to ISAJET (for msugra. x only) and/or 
SOFTSUSY within respectively the AMSB, GMSB, mSUGRA or NUHM parameter spaces. 

test_modelef f and test_standmod calculate the relic density, using the parameters con- 
tained in the SLHA file whose name has to be passed as input parameter, in the cosmological 
models described in the Appendices. 



4 Input and output description 

The input and output of the main programs are described in the following. 
4.1 mSUGRA inputs 

The program msugra . x computes the observables in the mSUGRA parameter space, using 
ISAJET and/or SOFTSUSY to generate the mass spectra. If only one of these generators is 
available the corresponding #define in msugra. c has to be commented. The necessary 
arguments to this program are: 

• mo: universal scalar mass at GUT scale, 

• mi/2- universal gaugino mass at GUT scale, 

• ^0- trilinear soft breaking parameter at GUT scale, 

• tan (3: ratio of the two Higgs vacuum expectation values. 

Optional arguments can also be given: 

• sign{fi): sign of Higgsino mass term, positive by default, 

• mf°'^: top quark pole mass, by default 172.4 GeV, 

• mb{rnf,): scale independent b-quark mass, by default 4.2 GeV (option only available 
for SOFTSUSY), 
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• as{Mz)- strong coupling constant at the Z-boson mass, by default 0.1176 (option 
only available for SOFTSUSY). 

If the arguments are not specified, a message will describe the needed parameters in the 
correct order. 

With SOFTSUSY 2.0.18 and ISAJET 7.78, running the program with: 

./msugra.x 1000 1000 10 

returns 

SLHA file generated by SOFTSUSY 
delta0=8. 116e-02 
BR_bsgamma=3 . 082e-04 
BR_Btaunu=l . lOOe-04 
Rtaunu=9.980e-01 
BR_Kniunu/BR_pimunu=6 . 454e-01 
R123=1.000e+00 
BR_BDtaunu=6 . 976e-03 
BR_BDtaunu/BR_BDenu=2 . 974e-01 
BR_Bsmumu=3 . 095e-09 
BR_Dstaunu=4 . 818e-02 
BR_Dsmunu=4 . 975e-03 
a_muon=9 . 756e-ll 
charged_LSP=0 
excluded_iiiass=0 
omega=9 . 087e+00 

SLHA file generated by ISAJET 
delta0=8.094e-02 
BR_bsgamma=3 . 072e-04 
BR_Btaunu=l . 099e-04 
Rtaunu=9.980e-01 
BR_Kiiiunu/BR_pimunu=6 . 454e-01 
R123=1.000e+00 
BR_BDtaunu=6 . 968e-03 
BR_BDtaunu/BR_BDenu=2 . 974e-01 
BR_Bsmumu=3 . 468e-09 
BR_Dstaunu=4 . 813e-02 
BR_Dsmunu=4 . 970e-03 
a_muon=1.013e-10 
charged_LSP=0 
excluded_iiiass=0 
omega=9. 102e+00 

where deltaO refers to the isospin symmetry breaking in S — > K*j decays, BR_bsgamma 
the branching ratio of — > Xs^, BR_Btaunu the branching ratio of Bu — > tv-j-j Rtaunu the 



12 



normalized ratio to the SM value, BR_Kmunu/BR_pimunu the ratio BR(ir /i;y^)/BR(7r — > 
fJ^i'n), R123 the ratio i?£23j BR_BDtaunu the branching ratio of B ^ D^tUt, BR_BDtaunu/BR_BDenu 
the ratio BR{B L>°ri/^)/BR(B eue), BR_Bsinumu the branching ratio of Bg — > , 

BR_Dstaunu and BR_Dsmunu the branching ratios of Ds rVj- and Ds fii^fi respectively, 
a_muon the deviation in the anomalous magnetic moment of the muon, charged_LSP in- 
dicates that the LSP is charged if equal to 1, excluded_mass that the point is already 
excluded by the direct searches if it is equal to 1, and omega is the relic density O/i^. More 
details on the definitions and calculations of the flavor observables are given in [11] . 

4.2 AMSB inputs 

The program amsb . x computes the observables using the corresponding parameters gener- 
ated by SDFTSUSY in the AMSB scenario. The necessary arguments to this program are: 

• ?Tio: universal scalar mass at GUT scale, 

• ^3/2- gravitino mass at GUT scale, 

• tan/3: ratio of the two Higgs vacuum expectation values. 

Optional arguments are the same as for mSUGRA. If the input parameters are missing, a 
message will ask for them. 

With SOFTSUSY 2.0.18, running the program with: 

./amsb.x 500 5000 5 -1 

returns 

delta0=7.669e-02 
BR_bsgainina=3 . 296e-04 
BR_Btaunu=l . 097e-04 
Rtaunu=9.953e-01 
BR_Kmunu/BR_pimunu=6 . 454e-01 
R123=1.000e+00 
BR_BDtaunu=6 . 971e-03 
BR_BDtaunu/BR_BDenu=2 . 972e-01 
BR_Bsmumu=3 . 349e-09 
BR_Dstaunu=4 . 818e-02 
BR_Dsmunu=4 . 975e-03 
a_muon=-6 . 469e-10 
excluded_iiiass=l 
oiiiega=1.291e-04 

4.3 GMSB inputs 

The program gmsb . x computes the observables using the GMSB parameters generated by 
SOFTSUSY. The necessary arguments to this program are: 

• A: scale of the SUSY breaking in GeV (usually 10000-100000 GeV), 

• Mmess- messenger mass scale (> A), 
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• N^: equivalent number of 5 + 5 messenger fields, 

• tan/3: ratio of the two Higgs vacuum expectation values. 

Optional arguments are the same as for mSUGRA, with an additional one: 

• CGrav 1): ratio of the gravitino mass to its value for a breaking scale A, 1 by 
default. 

Again, in the case of lack of arguments, a message will be displayed. 

With SOFTSUSY 2.0.18, running the program with: 
./gmsb.x 2e4 5e6 1 10 

returns 

delta0=6.309e-02 
BR_bsgamnia=4 . 285e-04 
BR_Btaunu=8 . 078e-05 
Rtaunu=7.331e-01 
BR_Kniunu/BR_piiiiunu=6 . 439e-01 
R123=9.988e-01 
BR_BDtaunu=6 . 542e-03 
BR_BDtaunu/BR_BDenu=2 . 789e-01 
BR_Bsmumu=3 . 802e-09 
BR_Dstaunu=4 . 806e-02 
BR_Dsmunu=4 . 962e-03 
a_muon=3 . 778e-08 
excluded_mass=l 
omega=3 . 497e-02 

4.4 NUHM inputs 

The program nuhm . x computes the observables using the NUHM parameters generated by 
SOFTSUSY. The necessary arguments to this program are the same as for mSUGRA, with 
two additional ones, the values of /U and tua- 

• itiq: universal scalar mass at GUT scale, 

• 1711/2: universal gaugino mass at GUT scale, 

• Aq: trilinear soft breaking parameter at GUT scale, 

• tan/3: ratio of the two Higgs vacuum expectation values, 

• /i: /U parameter, 

• niA'- CP-odd Higgs mass. 

Optional arguments can also be given: 

• m^"^'^: top quark pole mass, by default 172.4 GeV, 

• rnh{rnb): scale independent b-quark mass, by default 4.2 GeV, 
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• as{Mz)- strong coupling constant at the Z-boson mass, by default 0.1176. 
In the absence of arguments, a message will be shown. 

With SOFTSUSY 2.0.18, running the program with: 
./nuhm.x 500 500 50 500 500 

returns 

deltaO=l .069e-01 
BR_bsgamma=2 . 006e-04 
BR_Btaunu=6 . 775e-05 
Rtaunu=6. 148e-01 
BR_Kmunu/BR_pimunu=6 . 430e-01 
R123=9.982e-01 
BR_BDtaunu=6 . 333e-03 
BR_BDtaunu/BR_BDenu=2 . 700e-01 
BR_Bsmumu=3 . 338e-08 
BR_Dstaunu=4 . 796e-02 
BR_Dsiiiunu=4 . 952e-03 
a_niuon=2 . 188e-09 
excluded_mass=0 
oiiiega=7 . 533e-02 

4.5 SLHA input file 

The program slha.x calculates the observables while reading the needed parameters in a 
given SLHA file. For example, the command 

./slha.x example. lha 

returns 

delta0=8.259e-02 
BR_bsgamma=2 . 974e-04 
BR_Btaunu=l . 097e-04 
Rtaunu=9.966e-01 
BR_Kmunu/BR_piiiiunu=6 . 454e-01 
R123=1.000e+00 
BR_BDtaunu=6 . 966e-03 
BR_BDtaunu/BR_BDenu=2 . 973e-01 
BR_Bsmumu=3 . 470e-09 
BR_Dstaunu=4 . 813e-02 
BR_Dsmunu=4 . 969e-03 
a_muon=l .938e-10 
excluded_mass=0 
omega=l . 188e+01 

If the SLHA file provided to slha.x is invalid, a message will be displayed: 
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• Invalid point means that the SLHA generator had not succeeded in generating the 
mass spectrum {e.g. due to the presence of tachyonic particles). 

• Model not yet implemented means that the SLHA file is intended for a model not 
implemented in Superlso, such as i?-parity violating models. 

• Invalid SLHA file means that the SLHA file is broken and misses important pa- 
rameters. 



4.6 Alternative QCD equations of state 

The program test_modelef f . x calculates the relic density while reading the needed pa- 
rameters in the SLHA file, for the different QCD equations of state {i.e. alternative models 
of QcS and /icff) described in Appendix O For example, the command 

. /test_modelef f . X example. lha 

returns 

Dependence of the relic density on the calculation of heff and geff 

For model_eff=l (model A): omega=l . 188e+01 

For model_eff=2 (model B (default)): omega=l . 188e+01 

For model_eff=3 (model B2) : omega=l . 196e+01 

For model_eff=4 (model B3) : omega=l . 181e+01 

For model_eff=5 (model C) : omega=l . 189e+01 

For model_eff=0 (old model): omega=l . 165e+01 



4.7 Effective energy and entropy densities 

The program test_standmodel . x reads the needed parameters in the SLHA file, and calcu- 
lates the relic density while adding to the standard cosmological model an effective energy 
density such that 

PD = KpPrad{TBBN){T /TbEnY" , (2) 

and/or an effective entropy density 

SD = l^sSrad{TBBN){T /TbBnY' , (3) 

which modify the Early Universe properties without having observational consequences if 
chosen adequately ^Q\. A description of the model and of the related equations can be 
found in Appendix [Dj The necessary arguments to this program are^i 

• SLHA file name, 

• Kp-. ratio of dark energy density over radiation energy density at BBN time (prefer- 
entially < 1), 

• Hp-, dark energy density decrease exponent (preferentially > 4), 

• Kg ■ ratio of dark entropy density over radiation entropy density at BBN time (prefer- 
entially < 1), 

"'■The preferential values given inside the brackets correspond to cosmological models without observational 
consequences, i.e. as valid as the cosmological standard model. 
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Figure 1: Wes in function of pefT, computed with Superlso Relic (dashed green hne), and 
with DarkSusy (red crosses). This comparison shows an excellent agreement. 

• Hg'. dark entropy density decrease exponent (preferentially > 3). 

Note that Up = 4 corresponds to a radiation-like energy density, = 6 to a quintessence- 
like energy density and = 8 to a decaying scalar field energy density. Also, = 3 
corresponds to a radiation-like entropy density and rig = 4 can appear in reheating models. 
For example, the command 

./test_standmod.x example. lha le-3 6 le-3 4 

returns 

For the cosmological standard model: 
omega=l . 188e+01 

For the specified model with dark density/entropy : 
omega=4.633e+02 

Using the aforementioned main programs as examples, the user is encouraged to write 
his/her own programs in order to, for example, perform scans in a given supersymmetric 
scenario, or test different cosmological models. 



5 Results 

Superlso Relic computes the relic density, and the results have been compared extensively 
to those of DarkSusy and Micromegas. A very good agreement has been found even at the 
level of the calculation of the effective annihilation rate Weg (see Appendix [A]), as can be 
seen in Fig. [TJ In general, the results of DarkSusy, Micromegas, and Superlso Relic differ 
only by a few percents, but in some rare cases where a Higgs resonance occurs approximately 



DarkSusy 
Superlso Relic 
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Figure 2: Constraints on the NUHM parameter plane (hjUia), in the standard cosmological 
model (left), and in presence of a tiny energy overdensity with Kp = 10~^ and Up = 6 
(right). The red points are excluded by the isospin asymmetry of i? — > K*^, the gray area 
is excluded by direct collider limits, the yellow zone involves tachyonic particles, and the 
blue strips are favored by the WMAP constraints. 



at twice the mass of the LSP, the differences can be large. To avoid this problem, a very 
precise calculation of the widths of the Higgs bosons is required, and we decided to use the 
two-loop calculations of FeynHiggs to obtain a better evaluation of the relic density in this 
kind of scenarios. 

Superlso Relic can also be used in order to constrain SUSY parameter spaces, as it 
provides many different observables from flavor physics as well as the relic density. It allows 
in particular to test easily the influence of the cosmological model by modifying for example 
the QCD equation-of-state (Appendix [C]) or the expansion rate (Appendix [D]) , as can be 
seen in Fig. [2l This unusual feature will be further developed in the next versions of the 
program. 
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Appendix A Thermally averaged annihilation cross section 



The computation of the thermally averaged annihilation cross section {av) is the most time 
consuming part of the relic density computation, as it requires the computation of the 
many annihilation and co-annihilation amplitudes. One can define the annihilation rate of 
supersymmetric particles i and j into SM particles k and I |18^ I19j: 



ij^kl 



Pkl 



IQTi'^gigjSkiy/s 



E 



internal d.o.f. 



\M{ij kl)\^dn , 



(4) 



where is the transition amplitude, s is the center-of-mass energy, gi is the number of 
degrees of freedom of the particle i, pki is the final center-of-mass momentum such as 



[s - {ruk + mi)"^] [s - {ruk - mif] 
27i 



Pkl 



211/2 



(5) 



Ski is a symmetry factor equal to 2 for identical final particles and to 1 otherwise, and the 
integration is over the outgoing directions of one of the final particles. Moreover, an average 
over initial internal degrees of freedom is performed. 



We also define an effective annihilation rate VFcg by 



(6) 



with 

Practically, we compute 
dcos 9 



2 - Amlgp . 



E 

ijkl 



PijPkl 



^TTgLSpPcsSkiVs 



E 

helicities 



diagrams 



kl) 



(7) 



(8) 



where 9 is the angle between particles i and k. We integrate over cos 9 numerically by means 
of the Gauss-Legendre quadrature of order 5. 

Since Wefr(\/s) does not depend on the temperature T, it can be tabulated once for each 
model point. It is however important to make sure that the maximum y/s in the table is 
large enough to include all important resonances, thresholds and coannihilation thresholds. 

To improve the calculation speed, we use two different thresholds: 

• a cut such that the coannihilation of SUSY particles i and j is only taken into account 
if 

rtli + rrij < ^/s^ut coann ) (9) 

where we have taken 



^cut coann 



3 misp 



(10) 



19 



• a maximum energy up to which Wes{^/s) is calculated, such that 



'^niLSP - Tfo log(Se) 



ill) 



where Tf^ = 25 GeV is a typical upper limit freeze-out temperature, and is the 
Boltzmann suppression factor limit that we fixed at 10~^ j20] . 

The thermal average of the effective cross section is then 



mlspT 



E 



Qi m- /rrii 

2-^2 

QLSP ml 



2 ' 



(12) 



where Ki and K2 are the modified Bessel functions of the second kind of order 1 and 2 
respectively. The average is performed numerically using a Gaussian integration, and the 
CO limit can be safely replaced by Pcs[\fSrnax) using the properties of Ki. 



Appendix B Cosmological Standard Model 

The cosmological standard model is based on a Priedmann-Lemaitre Universe filled with 
radiation, baryonic matter and cold dark matter, approximately flat and incorporating 
a cosmological constant accelerating its expansion. Before recombination, the Universe 
expansion was dominated by a radiation density, and therefore the expansion rate H of the 
Universe is determined by the Friedmann equation 

= -^Prad , (13) 

where ^ 

PradiT) = g^^m^T" (14) 

is the radiation density and g^^^ is the effective number of degrees of freedom of radiation. 
The computation of the relic density is based on the solution of the Boltzmann evolution 
equation [151 [T9] 

dn/dt = —3Hn — {a^Qv){n'^ — n^q) , (15) 

where n is the number density of all supersymmetric particles, neq their equilibrium den- 
sity, and (cTggv) is the thermal average of the annihilation rate of the supersymmetric 
particles to the Standard Model particles. By solving this equation, the density number of 
supersymmetric particles in the present Universe and consequently the relic density can be 
determined. 

The ratio of the number density to the radiation entropy density, Y{T) = n{T) / s{T) can 
be defined, where 

s{T) = K^{T)—T\ (16) 
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/igg- is the effective number of entropic degrees of freedom of radiation. Combining Eqs. 
()13p and (|15p and defining x = m^gp/T, the ratio of the LSP mass over temperature, yield 



with 

The freeze-out temperature Tf is the temperature at which the LSP leaves the initial thermal 
equilibrium when Y{Tf) = (1 + 5)YQq(Tf), with 5 ~ 1.5. The relic density is obtained by 
integrating Eq. (fTTj) from x = to mj^gp/To, where Tq = 2.726 K is the temperature of 
the Universe today [IBl [l9] : 

^ m^sMTomTo)h^ ^ 2.755 x lO'^Yin) , (19) 
where is the critical density of the Universe, such as 



Hi = ^-^pI , (20) 



Hq being the Hubble constant. 



In practice, to improve the speed of the code, the freeze-out temperature Tj is determined 
by solving the implicit equation: 

and the evolution equation (I15p is only solved from T = Tf to To, with the initial condition 
Y{Tf) = (1 + 5)YQq. This method is known to provide results with less than a few percent 
error for the calculation of the relic density. 



Appendix C QCD equation of state 

To evaluate the relic density, it is necessary to know the number of effective degrees of 
freedom g^Q and h^Q which give access to the energy and entropy densities of radiation. To 
compute them, one generally assumes that above the QCD phase transition critical tem- 
perature Tc ~ 200 MeV, the primordial plasma is weakly interacting because of asymptotic 
freedom, and can therefore be treated as an ideal gas. 

However, non-perturbative studies have shown that the QCD plasma departs strongly from 
the ideal gas behavior at high temperatures, and more realistic models have been studied 
in [T7] . In these models, below Tc the hadronic degrees of freedom are modeled by an inter- 
acting gas of hadrons, while above Tc the quarks and gluons are taken to interact and are 
replaced by hadronic models. In Superlso Relic, the models depicted in [T7] are available, 
and can be selected in the routine Init_modelef f (int niodel_ef f , struct relicparam* 
paramrelic) by setting the value of model_eff (see subsection 12. 2p . The different models 
are: 
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• Model A (model_ef f =1) which ignores hadrons completely. 

• Model B (model_ef f =2) which considers Tc = 154 MeV, and models hadrons as a 
gas of free mesons and hadrons, with a sharp switch to the hadronic gas at T^g = Tc- 

• Model B2 (model_ef f =3) is a variation of model B constructed by scaling the lattice 
data by 0.9. 

• Model B3 (niodel_ef f =4) is a variation of model B constructed by scaling the lattice 
data by 1.1. 

• Model C (model_ef f =5) which assumes Tc = 185.5 MeV, and models hadrons as a 
gas of free mesons and hadrons, with a sharp switch to the hadronic gas at T^g = 200 
MeV. 

• Old Model (model_ef f =0) is the old model with an ideal gas. 

An example main program is provided as test_modelef f . c. For more information about 
these models, the reader is referred to jl7j . 



Appendix D Modified Cosmological Model 



^ = -2,Hn - {av) (n^ - n%) + Nd , (22) 



The density number of supersymmetric particles is determined by the Boltzmann equation: 

dn 
Itt 

where n is the number density of supersymmetric particles, {av) is the thermally averaged 
annihilation cross-section, H is the Hubble parameter, rieq is the relic particle equilib- 
rium number density. The term No has been added to provide a parametrization of the 
non-thermal production of SUSY particles. The expansion rate H is determined by the 
Friedmann equation: 

H^=^-^{Prad + PD), (23) 

where prad is the radiation energy density, which is considered as dominant before BBN in 
the standard cosmological model. Following [9l|T0], /O^ is introduced as an effective dark 
density which parametrizes the expansion rate modification. The entropy evolution reads: 

^ = -3Fs + Sb, (24) 
at 

where s is the total entropy density. S o parametrizes here effective entropy fluctuations due 
to unknown properties of the Early Universe. The radiation energy and entropy densities 
can be written as usual: 

y^2 27]-2 
Prad = 9cs{T)—T^ , Srad = Ks{T)-—T^ . (25) 

30 45 
Separating the radiation entropy density from the total entropy density, i.e. setting s = 
Srad + sd where sd is an effective entropy density, the following relation between sd and 
TiD can be derived: 



47r3G 



1 /icff ^C^SZ) 



3 „i/2 dT 



(26) 
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Following the standard relic density calculation method [THl [12], we introduce Y = n/s, 
and Eq. (j22p becomes 



dY _ rriLSP I tt 1/2 ( 1 + sd 



'/leff(T)^r3)' (1 + 5^)2 



45 

(27) 

where x = uilsp/T, ttilsp being the mass of the relic particle, and 

= ——-^^ , PD = , (28) 

and 

45 , 1 ■s-^ o^, /m,- 



y = y g .^,2^, ( ^ ) , (29) 

i 

where i runs over all supersymmetric particles of mass rrii and with gi degrees of freedom. 
Following the methods described in Appendix [Bl the relic density can then be calculated: 

Qh"^ = 2.755 X W^YoniLSp/GeV . (30) 

where Yq is the present value of Y. In the limit where pn = sd = ^d = ^D = 0, usual 
relations are retrieved. We should note here that sd and T,d are not independent variables. 

In Superlso Relic v2.5, we adopt the parametrizations described in [9l [10] for and 

SD- 

PD = l^pPrad{TBBN){T/TBBN)"''' (31) 

and 

Sd = l^sSrad{TBBN){T /Tbbn) \ (32) 

where Tbbn stands for the BBN temperature. K,p and Kg are respectively the ratio of 
effective dark energy/entropy density over radiation energy /entropy density, and Up and Ug 
are parameters describing the behavior of the densities. We refer the reader to [9l [10] for 
detailed descriptions and discussions of these parametrizations. 
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